Add buffers around points in Python

Adding buffers around points is essential for species distribution models. This post shows how to add buffers with different radii around the red and green Kangaroo’s Paw using Python.

Eukaryota
Plantae
Charophyta
Commelinales
Haemodoraceae
Chordata
Anigozanthos
Buffers
Python
Authors

Amanda Buyan

Dax Kellie

Published

November 23, 2023

Author

Amanda Buyan
Dax Kellie

Date

23 November 2023

Buffers are used in many ecological spaces. Most notably, they are used in Species Distribution Models to extend the possible range of rare or threatened species to create a better model. They can also be used to capture occurrences or areas which may have been neglected by a user’s original data.

In this post, we show how to use {galah-python}, {geopandas}, {shapely} and {matplotlib} to add a buffer to occurrences of Kangaroo Paw, specifically the Red and Green Kangaroo Paw.

Choosing the right taxon

Kangaroo Paw is a common name for a number of species, represented by the genus Anigozanthos. They are perennial plants, native to the south-west of Western Australia. They are unique, bird-attracting flowers which open at the apex with six claw-like structures resembling kangaroo paws, hence their name!

To begin, we will look for the number of occurrences for the genus Anigozanthos. To do this, we will first start by loading galah-python.

import galah
galah.galah_config(email="your-email-here")

Now, we will use the atlas_counts() function to see a count of the number of occurrences of all species of Kangaroo’s Paw in the ALA.

galah.atlas_counts(taxa="Anigozanthos")
   totalRecords
0            84

It is surprising to note that there are only 84 occurrences coming up for all species in this genus! To ensure that we have the correct name, we will use search_taxa() to see if there is anything going on with this genus.

galah.search_taxa(taxa="Anigozanthos")
                            scientificName                                   taxonConceptID      rank  kingdom      phylum         order         family         genus   issues
0  Anigozanthos 'Kings Park Western Glory'  https://id.biodiversity.org.au/name/apni/226062  unranked  Plantae  Charophyta  Commelinales  Haemodoraceae  Anigozanthos  noIssue

If you scroll across to the rank column that is returned by this function, unfortunately, the search item that comes up for search_taxa() is unranked, which explains why we are getting very few records! What we can do now to refine our search is add a higher order of taxonomic classification to our query. By searching on the ALA species pages, we know that Anigozanthos is a genus belonging to the Haemodoraceae family. galah-python has an argument to search_taxa() called scientific_name, which allows the user to disambiguate their query.

Now, if we add the extra taxonomic information, we get thte following query:

galah.search_taxa(scientific_name={"family": ["Haemodoraceae"],"genus": ["Anigozanthos"],"scientificName": ["Anigozanthos"]})
  scientificName scientificNameAuthorship                                    taxonConceptID   rank  kingdom      phylum         order         family         genus         vernacularName   issues
0   Anigozanthos                  Labill.  https://id.biodiversity.org.au/node/apni/2890004  genus  Plantae  Charophyta  Commelinales  Haemodoraceae  Anigozanthos  Australian Sword Lily  noIssue

Here, if we look at the information, we can see that our scientificName matches the genus we want, and the rank is now listed as genus. Now, if we get the total record count using this disambiguation:

galah.atlas_counts(scientific_name={"family": ["Haemodoraceae"],"genus": ["Anigozanthos"],"scientificName": ["Anigozanthos"]})
   totalRecords
0          5494

Another, shorter, way to write this query is by using the scientificName and scientificNameAuthorship as the taxon name. To test that we do get the same number of records, we can simply run

galah.atlas_counts(taxa="Anigozanthos Labill.")
   totalRecords
0          5494

Downloading occurrence records

Now that we know the correct genus, we can start to download occurrence records. Since there are ~6000 records, let’s choose a smaller subset to draw buffers around. Since a breeding program was started in 2007 by Kings Park Botanic and Garden Board to protect the Kangaroo Paw from disease and the impact of climate changes, let’s only include occurrence records starting from 2007 onwards.

galah.atlas_counts(taxa="Anigozanthos Labill.",filters="year>=2007")
   totalRecords
0          1930

1936 is a more manageable number. However, there are different ways to observe a species. Since we want to only include ones humans have observed in the wild, we will add the filter basisOfRecord=HUMAN_OBSERVATION

galah.atlas_counts(taxa="Anigozanthos Labill.",filters=["year>=2007","basisOfRecord=HUMAN_OBSERVATION"])
   totalRecords
0          1667

Now, we can download all occurrences of all species of Kangaroo Paw and plot it on a map to do an initial check that we have all records in Kangaroo Paw’s natural habitat, southwest WA.

import geopandas as gpd
import matplotlib.pyplot as plt
kanga_paw_occ = galah.atlas_occurrences(taxa="Anigozanthos Labill.",filters=["year>=2007","basisOfRecord=HUMAN_OBSERVATION"])
states = gpd.read_file("STE_2021_AUST_GDA2020.shp")
states = states.to_crs(4326)
ax = states.plot(edgecolor = "#5A5A5A", linewidth = 0.5, facecolor = "white", figsize = (12,6))
plt.scatter(kanga_paw_occ['decimalLongitude'],kanga_paw_occ['decimalLatitude'], c = "red")

Here, we can see that the ALA has records that are outside WA, and if we want to only draw a buffer region around species in Western Australia, we will have to filter out the points in the east. To do this, we will provide the polygon representing WA to atlas_occurrences() to remove the occurrences seen in the Eastern parts of Australia.

kanga_paw_occ_pol = galah.atlas_occurrences(
    taxa="Anigozanthos Labill.",
    filters=["year>=2007","basisOfRecord=HUMAN_OBSERVATION"],
    polygon = states[states["STE_NAME21"] == "Western Australia"]["geometry"][4],
    simplify_polygon=True,
)
ax = states[states["STE_NAME21"] == "Western Australia"].plot(edgecolor = "#5A5A5A", linewidth = 0.5, facecolor = "white", figsize = (12,6))
plt.scatter(kanga_paw_occ_pol['decimalLongitude'],kanga_paw_occ_pol['decimalLatitude'], c = "red")

Now, we have the occurrences we want.

galah.atlas_counts(
    taxa="Anigozanthos Labill.",
    filters=["year>=2007","basisOfRecord=HUMAN_OBSERVATION"],
    polygon = states[states["STE_NAME21"] == "Western Australia"]["geometry"][4],
    simplify_polygon=True,
    group_by="species",
    expand=False
)
                      species  count
0        Anigozanthos bicolor     67
1       Anigozanthos flavidus    117
2      Anigozanthos gabrielae      2
3        Anigozanthos humilis    456
4  Anigozanthos kalbarriensis      3
5      Anigozanthos manglesii    460
6       Anigozanthos preissii     14
7   Anigozanthos pulcherrimus      7
8          Anigozanthos rufus     92
9        Anigozanthos viridis     29

We can see that the most recorded species is Anigozanthos manglesii, which is the Red and Green Kangaroo Paw and the state emblem of WA. Since there are a decent number of records, but not as many as the entire genus, let’s choose to focus on the Red and Green Kangaroo Paw.

anigozanthos_manglesii = galah.atlas_occurrences(
    taxa = "Anigozanthos manglesii",
    filters=["year>=2007","basisOfRecord=HUMAN_OBSERVATION"],
    polygon = states[states["STE_NAME21"] == "Western Australia"]["geometry"][4],
    simplify_polygon=True
)

We can visualise these points to check that they are in the southwest corner of WA.

ax = states[states["STE_NAME21"] == "Western Australia"].plot(edgecolor = "#5A5A5A", linewidth = 0.5, facecolor = "white", figsize = (12,6))
plt.scatter(anigozanthos_manglesii['decimalLongitude'],anigozanthos_manglesii['decimalLatitude'], c = "red")

Adding Buffers

Now that we’ve got our occurrences, we can start adding buffers around these points. To do this, we will need the {shapely} package in Python. First, we will convert all of the decimalLongitude and decimalLatitude points to Point objects from {shapely}.

import shapely
from shapely.geometry import Point,Polygon
points_angiozanthos_manglesii = [Point((row['decimalLongitude'],row['decimalLatitude'])) for i,row in anigozanthos_manglesii[['decimalLongitude','decimalLatitude']].iterrows()]

Now that we have a list of all the points, we can add them to a geopandas DataFrame. This is so we can manipulate spatial data more easily than in pandas Dataframes. We will also set the Coordinate Reference System (CRS) to EPSG:4326, as this is the CRS of all ALA data points.

gdf_anigozanthos_manglesii = gpd.GeoDataFrame(anigozanthos_manglesii,geometry=points_angiozanthos_manglesii)
gdf_anigozanthos_manglesii.set_crs(epsg=4326, inplace=True)
     decimalLatitude  decimalLongitude             eventDate                          scientificName                                    taxonConceptID                              recordID                   dataResourceName occurrenceStatus                     geometry
0         -34.507979        116.907720  2021-10-30T06:55:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  57581a54-c592-4547-8864-f0f2034c91cf              iNaturalist Australia          PRESENT  POINT (116.90772 -34.50798)
1         -34.437661        116.954487  2022-10-29T05:53:16Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  4097f747-ab3d-44da-bba6-98bedaf094ba              iNaturalist Australia          PRESENT  POINT (116.95449 -34.43766)
2         -34.392112        116.033334  2018-10-14T12:17:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  d6309f9d-b324-4df2-b1b0-0d1d095ceba0              iNaturalist Australia          PRESENT  POINT (116.03333 -34.39211)
3         -34.383374        116.093064  2018-10-24T09:45:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  014cae6e-f46d-4a99-98aa-86089d7c6172              iNaturalist Australia          PRESENT  POINT (116.09306 -34.38337)
4         -34.382520        116.099223  2018-09-08T09:31:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  b1ff2a41-fdfa-4a82-a341-58dbc44dacd8              iNaturalist Australia          PRESENT  POINT (116.09922 -34.38252)
..               ...               ...                   ...                                     ...                                               ...                                   ...                                ...              ...                          ...
455       -27.700000        114.200000  2021-06-06T17:12:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  fcb60604-f38a-40ca-a6e6-7ad2d6154ee1        Earth Guardians Weekly Feed          PRESENT  POINT (114.20000 -27.70000)
456       -27.679641        114.239427  2023-09-04T03:27:38Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  e875313b-83b4-4207-af36-f7a2d89b69f2              iNaturalist Australia          PRESENT  POINT (114.23943 -27.67964)
457       -27.607720        114.420170  2011-09-20T11:44:00Z  Anigozanthos manglesii subsp. quadrans  https://id.biodiversity.org.au/node/apni/2894068  de617c68-bb47-4403-b355-9974bfba6df0              iNaturalist Australia          PRESENT  POINT (114.42017 -27.60772)
458       -27.564815        114.425332  2016-09-10T04:27:00Z  Anigozanthos manglesii subsp. quadrans  https://id.biodiversity.org.au/node/apni/2894068  f9c14a93-baf9-41c1-b3a1-fe7a9709ffa7              iNaturalist Australia          PRESENT  POINT (114.42533 -27.56481)
459       -27.260710        114.049158  2022-08-20T16:00:00Z                  Anigozanthos manglesii  https://id.biodiversity.org.au/node/apni/2900921  2788b724-4bbe-4979-aef7-9a9b34617ae7  ALA species sightings and OzAtlas          PRESENT  POINT (114.04916 -27.26071)

[460 rows x 9 columns]

Now, to add buffers of a certain radius, we have to convert our current CRS (which represents coordinates in degrees) to a CRS that represents coordinates in meters. This is so so we can directly add buffers in meters. For our example, we will be using EPSG:3577, which is Australian Albers, and a widely used CRS in Australia when one needs to use meters as a unit. This is likely different around the world, so be sure to check what CRS is right for your area.

gdf_anigozanthos_manglesii_meters = gdf_anigozanthos_manglesii.to_crs(3577)

Now, we will be creating five different buffers around our chosen points: 5km, 10km, 15km, 20km, and 25km. For each buffer radius, we will be adding a circle of the chosen radius around each point. We will then put them in a geodataframe so we can easily convert the buffers back from meters into degrees so they agree with our original points. Then, we will perform something called a unary_union, which is a function in shapely that allows you to unify many shapes into one shape object. We will then be storing this in a dictionary.

buffer_shapes = {}
buffer_lengths = {"5km": 5000, "10km": 10000,"15km": 15000,"20km": 20000,"25km": 20000}
for length in buffer_lengths:
  buffers = [row["geometry"].buffer(buffer_lengths[length]) for i,row in gdf_anigozanthos_manglesii_meters.iterrows()]
  gdf_buffers = gpd.GeoSeries(buffers).set_crs(3577)
  gdf_buffers_degrees = gdf_buffers.to_crs(4326)
  union_buffers_degrees = shapely.unary_union(gdf_buffers_degrees)
  buffer_shapes[length] = union_buffers_degrees

Now, we can finally plot the buffers! Below is a code using the dictionary we just created to run a for loop over all the buffers to easily plot them, and also colour them in different colours.

ax = states[states["STE_NAME21"] == "Western Australia"].plot(edgecolor = "#5A5A5A", linewidth = 0.5, facecolor = "white", figsize = (12,6))
ax.set_ylim([-34,-27])
ax.set_xlim([114,119])
colors = ["red","orange","green","purple","black"]
plt.scatter(anigozanthos_manglesii['decimalLongitude'],anigozanthos_manglesii['decimalLatitude'], c = "blue", s = 2)
for i,length in enumerate(buffer_lengths):
  for j,geom in enumerate(buffer_shapes[length].geoms):
    if j==0:
      plt.plot(*geom.exterior.xy,c=colors[i],lw=0.5,label=length)
    else:
      plt.plot(*geom.exterior.xy,c=colors[i],lw=0.5)
plt.legend(fontsize=16)

Final thoughts

We hope this point has helped make the steps of adding buffers around points on a map clearer and easier. For more advanced mapping in Python, check out our ALA Labs article on how to map invasive species.

Expand for session info
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.0 (2023-04-21)
 os       macOS 14.2
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_AU.UTF-8
 tz       Australia/Sydney
 date     2023-12-14
 pandoc   3.1.1 @ /Applications/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 htmltools   * 0.5.5   2023-03-23 [1] CRAN (R 4.3.0)
 sessioninfo * 1.2.2   2021-12-06 [1] CRAN (R 4.3.0)

 [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /Users/buy003/Library/r-miniconda-arm64/envs/r-reticulate/bin/python
 libpython:      /Users/buy003/Library/r-miniconda-arm64/envs/r-reticulate/lib/libpython3.9.dylib
 pythonhome:     /Users/buy003/Library/r-miniconda-arm64/envs/r-reticulate:/Users/buy003/Library/r-miniconda-arm64/envs/r-reticulate
 version:        3.9.16 | packaged by conda-forge | (main, Feb  1 2023, 21:38:11)  [Clang 14.0.6 ]
 numpy:          /Users/buy003/Library/r-miniconda-arm64/envs/r-reticulate/lib/python3.9/site-packages/numpy
 numpy_version:  1.25.2

──────────────────────────────────────────────────────────────────────────────